########################################################################################################
#### REPLICATION SCRIPT: I'M COMING OUT! HOW VOTER DISCRIMINATION PRODUCES EFFECTIVE LGBTQ LAWMAKERS 
### Jacob M. Lollis, UVA & Mackenzie R. Dobson, UVA & University of Notre Dame
########################################################################################################
# Required packages for replication 

library(foreign)
library(tidyverse)
library(readstata13)
library(ggplot2)
library(haven)
library(scales)
library(readxl)
library(usmap)
library(patchwork)
library(dplyr)
library(tidyr)
library(maps)

########################################################################################################
# IN-TEXT FIGURE REPLICATION
########################################################################################################

### Figure 1: LGBTQ Legislators Are More Effective Lawmakers ### 

# Regression coefficients and their standard errors were exported from STATA as a CSV file.
# This file is accessible on DataVerse as "Lollis&Dobson.2024.Figure1.xlsx."
# The model estimated corresponds with Table 4.1 in the manuscript's appendix and replication .do file. 

## MUST USE FILE  "Lollis&Dobson.2024.Figure1.xlsx"

# Import coefficients 
lgbtq.sles.df = read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.Figure1.xlsx")
lgbtq.sles.data.1 = lgbtq.sles.df[,2:3]

# Dropping "b" and "se" labels from coefficient export from STATA
subset.lgbtq.sles = lgbtq.sles.data.1[-c(1), ]

# Create a vector of an indicator variable that labels the predictors 
preds = rep(c("LGBTQ", "Percent LGBTQ", "Female", "Black","Hispanic", "Race (other)", 
              "White", "Democrat", "Seniority", "Committee Chair", "In Majority Party",
              "In Governor's Party", "In Majority Leadership", "In Minority Leadership", "Distance from Median",
              "Speaker/President", "Term Limits", "Professionalism", "Vote Share"), len=19)

# Create a vector of an indicator variable that labels the DV
dv= rep(c("SLES"), len=19)

# Create a vector of an indicator variable that labels the coefficient of interest as a different color 
group = rep(c("primary", "control", "control", "control", "control", "control", "control", "control",
              "control", "control", "control", "control", "control", "control", "control", "control",
              "control", "control", "control"), each=1)

model.label = cbind.data.frame(subset.lgbtq.sles, preds, group, dv)

# Rename column labels
colnames(model.label)[1] = "estimate"
colnames(model.label)[2] = "se"
colnames(model.label)[3] = "predictor"
colnames(model.label)[4] = "group"
colnames(model.label)[5] = "dv"

# Format numeric columns to avoid scientific notation
model.label$estimate <- format(model.label$estimate, scientific = FALSE)
model.label$se <- format(model.label$se, scientific = FALSE)

# Plot 
model.label$predictor= factor(model.label$predictor ,
                              levels = c("Vote Share",
                                         "Professionalism",
                                         "Term Limits",
                                         "Speaker/President",
                                         "Distance from Median",
                                         "In Minority Leadership",
                                         "In Majority Leadership",
                                         "In Governor's Party",
                                         "In Majority Party",
                                         "Committee Chair",
                                         "Seniority",
                                         "Democrat",
                                         "White",
                                         "Race (other)",
                                         "Hispanic",
                                         "Black",
                                         "Female",
                                         "Percent LGBTQ",
                                         "LGBTQ"))

estimate.num.lqbtq.sles = as.numeric(model.label$estimate)
std.num.time.lqbtq.sles = as.numeric(model.label$se)

ymin.lgbtq.sles = estimate.num.lqbtq.sles - 1.96 *std.num.time.lqbtq.sles
ymax.lgbtq.sles = estimate.num.lqbtq.sles + 1.96 *std.num.time.lqbtq.sles


lgbtq.sles.plot = ggplot(model.label) +
  theme_bw() +
  geom_pointrange(
    aes(
      x = predictor,
      y = estimate.num.lqbtq.sles,
      ymin = ymin.lgbtq.sles,
      ymax = ymax.lgbtq.sles,
      color = group
    ),
    position = position_dodge(width=.6),
    fatten = 4
  ) +
  geom_hline(yintercept=0, color="gray60", linetype=2) +
  coord_flip() +
  theme(axis.title.y=element_blank()) +
  theme(strip.text.y=element_text(angle=0)) +
  theme(text=element_text(size=13)) +
  ylab('Coefficient estimate') +
  theme(axis.text.y=element_text(hjust=0)) +
  theme(legend.position="none") +  
  scale_color_manual(values = c("grey70", "black"), guide = guide_legend(reverse=F)) 

lgbtq.sles.plot


### Figure 2: Out Legislators Are More Effective Lawmakers ### 

# Regression coefficients and their standard errors were exported from STATA as a CSV file.
# This file is accessible on DataVerse as "Lollis&Dobson.2024.Figure2.xlsx."
# The model estimated corresponds with Table 4.2 in the manuscript's appendix and replication .do file. 

## MUST USE FILE  "Lollis&Dobson.2024.Figure2.xlsx"

# Import coefficients 
out.sles.df = read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.Figure2.xlsx")
out.sles.data.1 = out.sles.df[,2:3]

# Dropping "b" and "se" labels from coefficient export from STATA
subset.out.sles = out.sles.data.1[-c(1), ]

# Create a vector of an indicator variable that labels the predictors 
preds = rep(c("Out During Election", "Percent Out", "Female", "Black","Hispanic", 
              "White", "Democrat", "Seniority", "Committee Chair", "In Majority Party",
              "In Governor's Party", "In Majority Leadership", "In Minority Leadership", "Distance from Median",
              "Speaker/President", "Term Limits", "Professionalism", "Vote Share"), len=18)

# Create a vector of an indicator variable that labels the DV
dv= rep(c("SLES"), len=18)


# Create a vector of an indicator variable that labels the coefficient of interest as a different color 
group = rep(c("primary", "control", "control", "control", "control", "control", "control", "control",
              "control", "control", "control", "control", "control", "control", "control", "control",
              "control", "control"), each=1)

model.label = cbind.data.frame(subset.out.sles, preds, group, dv)


# Rename column labels
colnames(model.label)[1] = "estimate"
colnames(model.label)[2] = "se"
colnames(model.label)[3] = "predictor"
colnames(model.label)[4] = "group"
colnames(model.label)[5] = "dv"

# Format numeric columns to avoid scientific notation
model.label$estimate <- format(model.label$estimate, scientific = FALSE)
model.label$se <- format(model.label$se, scientific = FALSE)



# Plot 
model.label$predictor= factor(model.label$predictor ,
                              levels = c("Vote Share",
                                         "Professionalism",
                                         "Term Limits",
                                         "Speaker/President",
                                         "Distance from Median",
                                         "In Minority Leadership",
                                         "In Majority Leadership",
                                         "In Governor's Party",
                                         "In Majority Party",
                                         "Committee Chair",
                                         "Seniority",
                                         "Democrat",
                                         "White",
                                         "Hispanic",
                                         "Black",
                                         "Female",
                                         "Percent Out",
                                         "Out During Election"))


estimate.num.out.sles = as.numeric(model.label$estimate)
std.num.time.out.sles = as.numeric(model.label$se)

ymin.out.sles = estimate.num.out.sles - 1.96 *std.num.time.out.sles
ymax.out.sles = estimate.num.out.sles + 1.96 *std.num.time.out.sles


out.sles.plot = ggplot(model.label) +
  theme_bw() +
  geom_pointrange(
    aes(
      x = predictor,
      y = estimate.num.out.sles,
      ymin = ymin.out.sles,
      ymax = ymax.out.sles,
      color = group
    ),
    position = position_dodge(width=.6),
    fatten = 4
  ) +
  geom_hline(yintercept=0, color="gray60", linetype=2) +
  coord_flip() +
  theme(axis.title.y=element_blank()) +
  theme(strip.text.y=element_text(angle=0)) +
  theme(text=element_text(size=13)) +
  ylab('Coefficient estimate') +
  theme(axis.text.y=element_text(hjust=0)) +
  theme(legend.position="none") +  
  scale_color_manual(values = c("grey70", "black"), guide = guide_legend(reverse=F)) 

out.sles.plot


### Figure 3: LGBTQ Legislators from Red States Are More Effective Lawmakers ### 

# Regression coefficients and their standard errors were exported from STATA as a CSV file.
# This file is accessible on DataVerse as "Lollis&Dobson.2024.Figure3.xlsx."
# The model estimated corresponds with Table 4.3 in the manuscript's appendix and replication .do file. 

## MUST USE FILE  "Lollis&Dobson.2024.Figure3.xlsx"

# Import coefficients 
redstate.bluestate.df <- read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.Figure3.xlsx")
redstate.bluestate.data.1 = redstate.bluestate.df[,2:3]
redstate.bluestate.data.2 = redstate.bluestate.df[,4:5]

colnames(redstate.bluestate.data.2) = colnames(redstate.bluestate.data.1)

redstate.bluestate.data.3 = rbind.data.frame(redstate.bluestate.data.1, redstate.bluestate.data.2)

# Dropping "b" and "se" labels from coefficient export from STATA
subset.redstate.bluestate = redstate.bluestate.data.3[-c(1,21), ]

# Create a vector of an indicator variable that labels the predictors  
predictor.table = rep(c("LGBTQ", "Percent LGBTQ", "Female", "Black", "Hispanic",
                        "Race (other)", "White", "Democrat", "Seniority", "Committee Chair",
                        "In Majority Party", "In Governor's Party", "In Majority Leadership", "In Minority Leadership",
                        "Distance from Median", "Speaker/President", "Term Limits", "Professionalism", "Vote Share"), len=38)

# Create a vector of an indicator variable that labels the group
group.table = rep(c("Red State", "Blue State"), each=19)

# Create a vector of an indicator variable that labels the DV
dv.table = rep(c("SLES"), len=38)

model.label = cbind.data.frame(subset.redstate.bluestate, predictor.table, group.table, dv.table)

# Rename column labels
colnames(model.label)[1] = "estimate"
colnames(model.label)[2] = "se"
colnames(model.label)[3] = "predictor"
colnames(model.label)[4] = "group"
colnames(model.label)[5] = "dv"

# Plot 
redstate.bluestate.est = subset(model.label)

redstate.bluestate.est$predictor= factor(redstate.bluestate.est$predictor ,
                                         levels = c("Vote Share",
                                                    "Professionalism",
                                                    "Term Limits",
                                                    "Speaker/President",
                                                    "Distance from Median",
                                                    "In Minority Leadership",
                                                    "In Majority Leadership",
                                                    "In Governor's Party",
                                                    "In Majority Party",
                                                    "Committee Chair",
                                                    "Seniority",
                                                    "Democrat",
                                                    "White",
                                                    "Race (other)",
                                                    "Hispanic",
                                                    "Black",
                                                    "Female",
                                                    "Percent LGBTQ",
                                                    "LGBTQ"))


estimate.num.redstate.bluestate = as.numeric(redstate.bluestate.est$estimate)
std.num.redstate.bluestate = as.numeric(redstate.bluestate.est$se)

ymin.redstate.bluestate = estimate.num.redstate.bluestate - 1.96 *std.num.redstate.bluestate
ymax.redstate.bluestate= estimate.num.redstate.bluestate + 1.96 *std.num.redstate.bluestate


redstate.bluestate.plot = ggplot(redstate.bluestate.est) +
  theme_bw() +
  geom_pointrange(
    aes(
      x = predictor,
      y = estimate.num.redstate.bluestate,
      ymin = ymin.redstate.bluestate,
      ymax = ymax.redstate.bluestate,
      color = group
    ),
    position = position_dodge(width=.6),
    fatten = 4
  ) +
  geom_hline(yintercept=0, color="gray60", linetype=2) +
  coord_flip() +
  theme(axis.title.y=element_blank()) +
  theme(strip.text.y=element_text(angle=0)) +
  theme(text=element_text(size=13)) +
  ylab('Coefficient estimate') +
  theme(axis.text.y=element_text(hjust=0)) +
  theme(legend.position="bottom") +
  scale_color_manual(
    name = "",
    guide = guide_legend(reverse=T) ,values = c("mediumblue","red3")) 


redstate.bluestate.plot

########################################################################################################
# APPENDIX FIGURE REPLICATION
########################################################################################################

#### Section 5: Additional Figures ####

### Figure 5.1: Overtime Trends in LGBTQ Lawmakers Being Out for All of Their Legislative Tenure: 1991-2017 ###

## MUST USE FILE  "Lollis&Dobson.2024.Figure5.1.xlsx"

always.out.df = read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.Figure5.1.xlsx")


always.out.plot = ggplot(data = always.out.df, aes(x = year, y = percent))+
  geom_line(color = "gray40", linewidth = 2) +
  theme_minimal() +
  scale_x_continuous(limits = c(1991, 2017), breaks = c(1991, 1993, 1995, 1997, 1999, 2001,
                                                        2003, 2005, 2007, 2009, 2011, 2013,
                                                        2015, 2017)) +
  labs(
    x = "Year",
    y = "Frequency (%) of LGBTQ Lawmakers Out") 


always.out.plot


### Figure 5.2: Overtime Trends in % LGBTQ and % Out Lawmakers Across State Legislatures: 1987-2018 ###

## MUST USE MULTIPLE .DTA FILES: 
# For 1980s-90s, use "Lollis&Dobson.2024.Figure5.2.pool1.dta"
# For 2000s, use "Lollis&Dobson.2024.Figure5.2.pool2.dta"
# For 2010s, use "Lollis&Dobson.2024.Figure5.2.pool3.dta"

# 80s & 90s
pool1 = read_dta("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.Figure5.2.pool1.dta")

# 2000s
pool2 = read_dta("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.Figure5.2.pool2.dta")

# 2010-2018
pool3 = read_dta("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.Figure5.2.pool3.dta")


# Aggregate percent LGBTQ and percent out from chamber level to legislature level 
agg.pool1 <- pool1 %>%
  group_by(state, term) %>%
  summarize(
    mean_pcLGBTQ = mean(pcLGBTQ),
    mean_pcOut = mean(pcOut),
  )

agg.pool2 <- pool2 %>%
  group_by(state, term) %>%
  summarize(
    mean_pcLGBTQ = mean(pcLGBTQ),
    mean_pcOut = mean(pcOut),
  )

agg.pool3 <- pool3 %>%
  group_by(state, term) %>%
  summarize(
    mean_pcLGBTQ = mean(pcLGBTQ),
    mean_pcOut = mean(pcOut),
  )

# 80s & 90s
map.8090.lgbtq = plot_usmap(regions = "state", data = agg.pool1, values = "mean_pcLGBTQ",  color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 4), breaks = seq(0, 4, by = 0.5),  # Keep increments at 0.5
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "% LGBTQ Lawmakers",
                        labels = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4),  # Specify the label format
                        expand = c(0, 0.5)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(0.55, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 10)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "1987-1999")


# Modify the linewidth of the borders
map.8090.lgbtq$layers[[2]]$aes_params$size <- 2

map.8090.lgbtq = map.8090.lgbtq +
  labs(title = "1987-1999") +
  theme(legend.position = 'right',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.5, "cm"),  # Increase the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

print(map.8090.lgbtq)



# 2000-2009
map.2000.lgbtq = plot_usmap(regions = "state", data = agg.pool2, values = "mean_pcLGBTQ", color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 7), breaks = seq(0, 7, by = 0.5),
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "",
                        labels = c(0, 0.5, 1, 1.5, seq(2, 7, by = 0.5)),
                        expand = c(0, 1)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(1, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 15)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "2000-2009")

# Modify the linewidth of the borders
map.2000.lgbtq$layers[[2]]$aes_params$size <- 2

map.2000.lgbtq = map.2000.lgbtq +
  labs(title = "2000-2009") +
  theme(legend.position = 'right',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.15, "cm"),  # Maintain the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),  # Increase the space between legend labels
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

print(map.2000.lgbtq)




# 2010 - 2018
map.1018.lgbtq = plot_usmap(regions = "state", data = agg.pool3, values = "mean_pcLGBTQ", color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 8), breaks = seq(0, 8, by = 0.5),
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "",
                        labels = c(0, 0.5, 1, 1.5, seq(2, 8, by = 0.5)),
                        expand = c(0, 1)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(1.15, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 15)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "2010-2018")

# Modify the linewidth of the borders
map.1018.lgbtq$layers[[2]]$aes_params$size <- 2

map.1018.lgbtq = map.1018.lgbtq +
  labs(title = "2010-2018") +
  theme(legend.position = 'bottom',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.5, "cm"),  # Maintain the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),  # Increase the space between legend labels
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

print(map.1018.lgbtq)

# Left column of the final figure
left = map.8090.lgbtq / map.2000.lgbtq / map.1018.lgbtq

left = left +  plot_layout(guides = "collect") & theme(legend.position = 'right') & theme(legend.text = element_text(size=11)) &
  theme(legend.title = element_text(size=17)) & theme(plot.title=element_text(size=17))

print(left)

# 80s & 90s
map.8090.lgbtq = plot_usmap(regions = "state", data = agg.pool1, values = "mean_pcOut",  color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 3.5), breaks = seq(0, 3.5, by = 0.5),  # Keep increments at 0.5
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "% Out-LGBTQ Lawmakers",
                        labels = c(0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5),  # Specify the label format
                        expand = c(0, 0.5)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(0.55, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 10)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "1987-1999")


# Modify the linewidth of the borders
map.8090.lgbtq$layers[[2]]$aes_params$size <- 2

map.8090.lgbtq = map.8090.lgbtq +
  labs(title = "1987-1999") +
  theme(legend.position = 'right',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.5, "cm"),  # Increase the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

print(map.8090.lgbtq)


# 2000-2009
map.2000.lgbtq = plot_usmap(regions = "state", data = agg.pool2, values = "mean_pcOut", color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 6), breaks = seq(0, 6, by = 0.5),
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "",
                        labels = c(0, 0.5, 1, 1.5, seq(2, 6, by = 0.5)),
                        expand = c(0, 1)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(1, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 15)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "2000-2009")

# Modify the linewidth of the borders
map.2000.lgbtq$layers[[2]]$aes_params$size <- 2

map.2000.lgbtq = map.2000.lgbtq +
  labs(title = "2000-2009") +
  theme(legend.position = 'right',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.15, "cm"),  # Maintain the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),  # Increase the space between legend labels
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

# Display the map
print(map.2000.lgbtq)


# 2010 - 2018
map.1018.lgbtq = plot_usmap(regions = "state", data = agg.pool3, values = "mean_pcOut", color = 'grey20', labels = TRUE) + 
  scale_fill_continuous(limits = c(0, 8), breaks = seq(0, 8, by = 0.5),
                        low = "grey70", high = 'grey10',
                        guide = "colorbar", na.value = "white", 
                        name = "",
                        labels = c(0, 0.5, 1, 1.5, seq(2, 8, by = 0.5)),
                        expand = c(0, 1)) +  # Adjust the spacing using expand
  geom_polygon(data = map_data("state"), aes(x = long, y = lat, group = group), color = "black", fill = NA, linewidth = 0.5) +  # Add state borders
  theme(legend.position = "right",
        legend.key.height = unit(1.15, "cm"),  # Adjust the legend key height
        legend.title = element_text(margin = margin(b = 15)),  # Adjust the margin between the title and the color bar
        panel.background = element_rect(colour = "black")) + 
  labs(title = "2010-2018")

# Modify the linewidth of the borders
map.1018.lgbtq$layers[[2]]$aes_params$size <- 2

map.1018.lgbtq = map.1018.lgbtq +
  labs(title = "2010-2018") +
  theme(legend.position = 'right',
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 14),  # Adjust the legend title size
        legend.key.width = unit(1.5, "cm"),  # Maintain the width of the color bar
        legend.spacing.x = unit(0.5, "cm"),  # Increase the space between legend labels
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5))  # Bigger and centered title

# Display the map
print(map.1018.lgbtq)


# Right column of the final figure



right = map.8090.lgbtq / map.2000.lgbtq / map.1018.lgbtq

right = right +  plot_layout(guides = "collect") & theme(legend.position = 'right') & theme(legend.text = element_text(size=11)) &
  theme(legend.title = element_text(size=17)) & theme(plot.title=element_text(size=17))
right


### Figure 5.3: Overtime Trends in LGBTQ Lawmakers by Identity Group: 1991-2017 ###

## MUST USE FILE  "Lollis&Dobson.2024.Figure5.3.xlsx"

id.overtime.df = read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.Figure5.3.xlsx")

id.overtime.df$id = factor(id.overtime.df$id, levels = c("Lesbian", "Gay", "Bisexual", "Trans", "Queer"))

id.overtime.plot = ggplot(data = id.overtime.df, aes(x = year, y = percent, group=id, color=id))+
  geom_line(linewidth=2) +
  scale_color_manual(values = c("red",
                                "darkorange",
                                "#FFCC00",
                                "green",
                                "blue")) +
  theme_minimal() +
  scale_x_continuous(limits = c(1991, 2017), breaks = c(1991, 1993, 1995, 1997, 1999, 2001,
                                                        2003, 2005, 2007, 2009, 2011, 2013,
                                                        2015, 2017)) +
  scale_y_continuous(limits=c(0, 1.2), breaks = c(0, 0.2, 0.4, 0.6, 0.8,
                                                  1, 1.2)) +
  labs(
       x = "Year",
       y = "Frequency (%) of Lawmakers by Identity Group") +
  guides(color = guide_legend(title = "Identity Group"))

id.overtime.plot

### Figure 5.4: Out LGBTQ Legislators Who Are Not First Time Seat Holders Are Less Effective Lawmakers ###

# Regression coefficients and their standard errors were exported from STATA as a CSV file.
# This file is accessible on DataVerse as "Lollis&Dobson.2024.Figure5.4.xlsx."
# The model estimated corresponds with Table 1 in the manuscript's main text and replication .do file. 

## MUST USE FILE  "Lollis&Dobson.2024.Figure5.4.xlsx"

# Import coefficients
replace.les.df <- read_xlsx("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson2024.Figure5.4.xlsx")
replace.les.data.1 = replace.les.df[,2:3]

# Dropping "b" and "se" labels from coefficient export from STATA
subset.replace.les = replace.les.data.1[-c(1), ]

# Create a vector of an indicator variable that labels the predictors  
predictor = rep(c("Not First LGBTQ from District", "Percent LGBTQ", "Female", "Black","Hispanic", 
                  "White", "Democrat", "Seniority", "Committee Chair", "In Majority Party",
                  "In Governor's Party", "In Majority Leadership", "In Minority Leadership", "Distance from Median",
                  "Speaker/President", "Term Limits", "Professionalism", "Vote Share"), len=18)


# Create a vector of an indicator variable that labels the group
group = rep(c("primary", "control", "control", "control", 
              "control", "control", "control", "control", 
              "control", "control", "control", "control",
              "control", "control", "control", "control",
              "control", "control"), len=18)

# Create a vector of an indicator variable that labels the DV
dv = rep(c("SLES"), len=18)

model.label = cbind.data.frame(subset.replace.les, predictor, group, dv)


# Rename column labels
colnames(model.label)[1] = "estimate"
colnames(model.label)[2] = "se"
colnames(model.label)[3] = "predictor"
colnames(model.label)[4] = "group"
colnames(model.label)[5] = "dv"

# Plot
replace.les.est = model.label

replace.les.est$predictor <- factor(replace.les.est$predictor,
                                    levels = c("Vote Share",
                                               "Professionalism",
                                               "Term Limits",
                                               "Speaker/President",
                                               "Distance from Median",
                                               "In Minority Leadership",
                                               "In Majority Leadership",
                                               "In Governor's Party",
                                               "In Majority Party",
                                               "Committee Chair",
                                               "Seniority",
                                               "Democrat",
                                               "White",
                                               "Hispanic",
                                               "Black",
                                               "Female",
                                               "Percent LGBTQ",
                                               "Not First LGBTQ from District"))

estimate.num.replace.les <- as.numeric(replace.les.est$estimate)
std.num.replace.les <- as.numeric(replace.les.est$se)

# Calculate ymin and ymax for 90% confidence intervals
ymin.replace.les <- estimate.num.replace.les - 1.645 * std.num.replace.les
ymax.replace.les <- estimate.num.replace.les + 1.645 * std.num.replace.les

replace.les.plot <- ggplot(replace.les.est) +
  theme_bw() +
  geom_pointrange(
    aes(
      x = predictor,
      y = estimate.num.replace.les,
      ymin = ymin.replace.les,
      ymax = ymax.replace.les,
      color = group
    ),
    position = position_dodge(width = .6),
    fatten = 3
  ) +
  geom_hline(yintercept = 0, color = "gray60", linetype = 2) +
  coord_flip() +
  theme(axis.title.y = element_blank()) +
  theme(strip.text.y = element_text(angle = 0)) +
  theme(text = element_text(size = 13)) +
  ylab('Coefficient estimate') +
  theme(axis.text.y = element_text(hjust = 0)) +
  theme(legend.position = "none") +
  scale_color_manual(
    values = c("primary" = "black", "control" = "gray60")
  )

replace.les.plot


#### Section 9: Transformation of the Dependent Variable ####

### Figure 9.1: Distribution of SLES and Z-Scored SLES ###

## MUST USE FILE  "Lollis&Dobson.2024.dta"

df = read_dta("/Users/mdobson2/Dropbox/LGBTQ Legislative Effectiveness/PS FINAL MANUSCRIPT/Data Replication/Lollis&Dobson.2024.dta")

# Original
hist.sles = ggplot(df, aes(x=SLES)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  theme_bw() +
  geom_rug(aes(SLES)) +
  geom_density(alpha=.2, fill="gray60") +
  ylab('Density') +
  xlab('Original SLES') + 
  scale_x_continuous(limits = c(-1,40))

hist.sles

# Transformed 
hist.z = ggplot(df, aes(x=SLES_z)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  theme_bw() +
  geom_rug(aes(SLES_z)) +
  geom_density(alpha=.2, fill="gray60") +
  ylab('Density') +
  xlab('Z-scored SLES') + 
  scale_x_continuous(limits = c(-1,40))

hist.z 


# Combine plots 
plot = hist.sles + hist.z
hist.plots = plot +  plot_layout(guides = "collect") & theme(legend.position = 'bottom')
hist.plots